Pump it Up Data Mining the Water Table
El objetivo es el de predecir la condición de un waterpoint para cada registro del conjunto de datos. Tenemos el siguiente conjunto de variables:
amount_tsh- cantidad de agua disponible para el waterpointdate_recorded- La fecha en que se ingresó el datofunder- Quién financió el pozogps_height- Altitud del pozoinstaller- Organización que instaló el pozolongitude- Coordenada GPS (longitud)latitude- Coordenada GPS (latitud)wpt_name- Nombre del waterpoint si existe algunonum_private-basin- Cuenca hidrográfica geográficasubvillage- Ubicación geográficaregion- Ubicación geográficaregion_code- Ubicación geográfica (codificada)district_code- Ubicación geográfica (codificada)lga- Ubicación geográficaward- Ubicación geográficapopulation- Población alrededor del pozopublic_meeting- Verdadero/Falsorecorded_by- Grupo que ingresa esta fila de datosscheme_management- Quién opera el punto de aguascheme_name- Nombre del esquema que opera el waterpointpermit- Si el punto de agua está permitidoconstruction_year- Año en que se construyó el waterpointextraction_type- Tipo de extracción que utiliza el waterpointextraction_type_group- Grupo de tipo de extracción que utiliza el waterpointextraction_type_class- Clase de tipo de extracción que utiliza el waterpointmanagement- Cómo se gestiona el waterpointmanagement_group- Grupo de gestión del waterpointpayment- Costo del aguapayment_type- Tipo de pago del aguawater_quality- Calidad del aguaquality_group- Grupo de calidad del aguaquantity- Cantidad de aguaquantity_group- Grupo de cantidad de aguasource- Fuente del aguasource_type- Tipo de fuente del aguasource_class- Clase de fuente del aguawaterpoint_type- Tipo de waterpointwaterpoint_type_group- Grupo de tipo de waterpoint
Los ficheros proporcionados son los siguientes:
train_values.csv: fichero con el conjunto de variables y observaciones con los que entrenar el modelo.train_labels.csv: fichero con las variables objetivo de cada una de las observaciones entrain_values.csv:- functional
- non functional
- functional needs repair
test_values.csv: fichero de prueba con el que maximizar las predicciones obtenidas con el correspondiente modelo
A lo largo del notebook, comentaré las distintas pruebas que he ido haciendo en el modelo y el cómo estas afectan al score de las predicciones del set de datos test.
0. Read the Data
train <- fread("./data/train.csv") %>% as.data.frame()
test <- fread("./data/test.csv") %>% as.data.frame()
trainlab <- fread("./data/train_labels.csv") %>% as.data.frame()
# put together train and test
TrainTest <- train %>% bind_rows(test) %>% as.data.table()1. Exploratory Data Analysis
paste0('Shape del dataframe: ',dim(train)[1],' filas y ',dim(train)[2],' columnas')## [1] "Shape del dataframe: 59400 filas y 40 columnas"
# Bar plot showing memory usage for each column
show_plot(inspect_mem(train))# Barplot of column types
show_plot(inspect_types(train))# Horizontal bar plot for categorical column composition
show_plot(inspect_cat(train))
Podemos ver que:
- Variables como
paymentypayment_typeoquantityyquantity_groupson idénticas, dado que tienen las mismas categorías y tienen el mismo volumen de datos por categoría. - La variable
recorded_bytan sólo tiene una categoría, lo cual no la hace interesante de cara al modelo. - Las variables
wpt_name,subvillage,ward,funder,installerodate_recordedtienen mucha cardinalidad - Algunas variables tienen una agrupación de categorías muy similares,
dado que pertenecen a un mismo grupo, como las
extraction.
# Bar plot of most frequent category for each categorical column
show_plot(inspect_imb(train))skim(train,amount_tsh,construction_year,district_code,gps_height,latitude,longitude,num_private,population,region_code)| Name | train |
| Number of rows | 59400 |
| Number of columns | 40 |
| _______________________ | |
| Column type frequency: | |
| numeric | 9 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| amount_tsh | 0 | 1 | 317.65 | 2997.57 | 0.00 | 0.00 | 0.00 | 20.00 | 350000.00 | ▇▁▁▁▁ |
| construction_year | 0 | 1 | 1300.65 | 951.62 | 0.00 | 0.00 | 1986.00 | 2004.00 | 2013.00 | ▅▁▁▁▇ |
| district_code | 0 | 1 | 5.63 | 9.63 | 0.00 | 2.00 | 3.00 | 5.00 | 80.00 | ▇▁▁▁▁ |
| gps_height | 0 | 1 | 668.30 | 693.12 | -90.00 | 0.00 | 369.00 | 1319.25 | 2770.00 | ▇▂▃▂▁ |
| latitude | 0 | 1 | -5.71 | 2.95 | -11.65 | -8.54 | -5.02 | -3.33 | 0.00 | ▃▅▅▇▂ |
| longitude | 0 | 1 | 34.08 | 6.57 | 0.00 | 33.09 | 34.91 | 37.18 | 40.35 | ▁▁▁▂▇ |
| num_private | 0 | 1 | 0.47 | 12.24 | 0.00 | 0.00 | 0.00 | 0.00 | 1776.00 | ▇▁▁▁▁ |
| population | 0 | 1 | 179.91 | 471.48 | 0.00 | 0.00 | 25.00 | 215.00 | 30500.00 | ▇▁▁▁▁ |
| region_code | 0 | 1 | 15.30 | 17.59 | 1.00 | 5.00 | 12.00 | 17.00 | 99.00 | ▇▁▁▁▁ |
# Histograms for numeric columns
show_plot(inspect_num(train))- La variable
amount_tshse distribuye en torno al0, aunque existen valores claramente atípicos que sugiere una posible transformación a tipo logarítmico o algún tipo de tratamiento de los datos. Lo mismo ocurre connum_privateypopulation. - La variable
construction_yeartiene un alto número de NAs, aunque no estén catalogados como tal. Podríamos imputarlos y crear una nueva variable que almacene si el dato era NA o no. Ocurre igual con la variabledistrict_code. - La variable
gps_heighttiene muchos valores en torno al 0. No sabemos si tiene sentido que las bombas están a una baja altura o no. - La variable
longitudetiene valores atípicos catalogados como 0, que posiblemente sean NAs
# Occurence of NAs in each column ranked in descending order
show_plot(inspect_na(train))# Correlation betwee numeric columns + confidence intervals
show_plot(inspect_cor(train))2. Preprocessing
2.1. Selecting categorical variables
# categorical variables df
cat <- TrainTest %>% select(where(is.character))
# get the cardinality of a category
levels <- data.frame("vars" = names(cat),
"levels" = apply(cat, 2,
function(x) length(unique(x))) )
# detele rownames
rownames(levels) <- NULL
# sort by cardinality
levels_sorted <- levels %>% arrange(levels)
levels_sorted %>% kbl() %>% kable_minimal()| vars | levels |
|---|---|
| recorded_by | 1 |
| source_class | 3 |
| management_group | 5 |
| quantity | 5 |
| quantity_group | 5 |
| quality_group | 6 |
| waterpoint_type_group | 6 |
| extraction_type_class | 7 |
| payment | 7 |
| payment_type | 7 |
| source_type | 7 |
| waterpoint_type | 7 |
| water_quality | 8 |
| basin | 9 |
| source | 10 |
| management | 12 |
| scheme_management | 13 |
| extraction_type_group | 13 |
| extraction_type | 18 |
| region | 21 |
| lga | 125 |
| ward | 2098 |
| funder | 2141 |
| installer | 2411 |
| scheme_name | 2869 |
| subvillage | 21426 |
| wpt_name | 45684 |
Variables con más de 3000 categorías pueden ser a priori
contraproducentes a la hora de crear un modelo, por lo que a priori
prescindiremos de ellas.
# select only categories with cardinality > 1 and < 2200
catvar <- levels %>%
filter(levels < 3000, levels > 1) %>%
select(vars)
# new dataframe
cat <- subset(cat, select = catvar$vars)Nos quedamos también con las variables con pocas categorías para posteriormente hacer dummies
# select only categories with cardinality > 1 and < 2200
dum <- levels %>%
filter(levels < 20, levels > 1) %>%
select(vars)
# new dataframe
dummies <- subset(cat, select = dum$vars)2.2. Selecting numerical variables
Una vez elegidas las variables categóricas, nos quedamos con las numéricas y las lógicas para hacer un nuevo dataset con dichas variables.
# numerical variables
num <- TrainTest %>% select(where(is.numeric))
# logical variables
logi <- TrainTest %>% select(where(is.logical))# join three dataframes
traintest <- cbind(num, logi,cat) %>% as.data.table()
traintest <- traintest %>% mutate(across(where(is.logical), ~ as.numeric(.x)))Una vez tenemos creado el dataset,fijémonos en algunas de las variables:
table(traintest$payment_type)##
## annually monthly never pay on failure other per bucket unknown
## 4570 10397 31712 4842 1314 11266 10149
table(traintest$payment)##
## never pay other pay annually
## 31712 1314 4570
## pay monthly pay per bucket pay when scheme fails
## 10397 11266 4842
## unknown
## 10149
# quantity_group: mismas categorias que con la variable quantity
table(traintest$quantity_group)##
## dry enough insufficient seasonal unknown
## 7782 41522 18896 5075 975
# quantity_group: mismas categorias que con la variable quantity
table(traintest$quantity)##
## dry enough insufficient seasonal unknown
## 7782 41522 18896 5075 975
Eliminamos por tanto las variables y los datos duplicados, ya que no
aportan información al modelo. Eliminamos también la variable
installer, ya que no conseguimos sacarle partido en el
modelo.
print_duplicates <- function(data) {
duplicated_rows <- duplicated(data)
duplicated_data <- data[duplicated_rows, ]
duplicated_data %<>% as.data.frame()
print(dim(duplicated_data)[1])
}
# duplicated rows
print_duplicates(traintest)## [1] 0
no tenemos valores duplicados
traintest <- traintest %>%
select(-quantity_group) %>%
select(-payment_type) %>%
select(-installer)2.3. Data Imputation
Como hemos comentado en la primera sección, las variables
construction_year, district_code,
longitude y gps_height tienen valores
missing catalogados como 0. En una primera
aproximación, hemos probado a pasar a NA y posteriormente
imputarlos,creando además una variable binaria que contabilice los datos
que eran NA:
# Function to impute zeros using missRanger method
impute <- function(data, columns, k, num_trees, impute = TRUE, seed = 0) {
for (col in columns) {
new_col <- paste0("NA_", col)
# create new variable that says if the value was NA or not
# if NA` -`> 1 ; else -> 0
data[[new_col]] <- ifelse(data[[col]] == 0, 1, 0)
data[[col]] <- ifelse(data[[col]] == 0, NA, data[[col]])
}
# now we impute values
if (impute == TRUE) {
data <- data %>%
# create variables that tell us if the row was NA or not
mutate(NA_public_meeting = ifelse(is.na(public_meeting),0,1)) %>%
mutate(NA_permit = ifelse(is.na(permit),0,1))
data_imp <- missRanger(data,
pmm.k = 5,
seed = 0,
maxiter = 100)
#data_imp <- kNN(data)#,k = floor(sqrt(dim(data)[1])))
}
return(data_imp)
}Pero la realidad es que el hecho de añadir las variables binarias al
dataset no mejora en absoluto el score del modelo. Las únicas
variables binarias que sí mejoran el rendimiento del modelo son las de
public_meeting y permit, por lo que optamos
por eliminar aquellas que empeoran el score. También he probado
distintos métodos de imputación de datos, como el MICE y
KNN:
# Function to impute zeros using KNN method
impute <- function(data, columns, k, num_trees, impute = TRUE, seed = 0) {
# now we impute values
if (impute == TRUE) {
data <- data %>%
# create variables that tell us if the row was NA or not
mutate(NA_public_meeting = ifelse(is.na(public_meeting),0,1)) %>%
mutate(NA_permit = ifelse(is.na(permit),0,1))
data_imp <- kNN(data)#,k = floor(sqrt(dim(data)[1])))
}
return(data_imp)
}# Function to impute zeros using MICE method
impute <- function(data, columns, k, num_trees, impute = TRUE, seed = 0) {
# now we impute values
if (impute == TRUE) {
data <- data %>%
# create variables that tell us if the row was NA or not
mutate(NA_public_meeting = ifelse(is.na(public_meeting),0,1)) %>%
mutate(NA_permit = ifelse(is.na(permit),0,1))
data_imp <- mice(data, m = 20, method = "pmm", seed = 0)
data_imp <- complete(data_imp)
}
return(data_imp)
}Pero ninguna de ellas mejora el score del modelo, por lo que optamos
por el método missRanger visto en clase.
# Funcion para la imputacion de ceros con missRanger
impute <- function(data, columns, k, num_trees, impute = TRUE, seed = 0) {
# now we impute values
if (impute == TRUE) {
data <- data %>%
# create variables that tell us if the row was NA or not
mutate(NA_public_meeting = ifelse(is.na(public_meeting),0,1)) %>%
mutate(NA_permit = ifelse(is.na(permit),0,1))
data_imp <- missRanger(data,
pmm.k = 5,
seed = 0,
maxiter = 100)
}
return(data_imp)
}
# impute data
traintestimp <- impute(traintest,
columns = c('construction_year', 'district_code', 'longitude','gps_height'),
k = 5,
num_trees = 100)2.4. Feature Engineering
Vamos con la sección más importante de todas, la creación de nuevas variables que sean capaces de mostrar información que anteriormente parecía oculta para el modelo. Para gran parte de la creación de variables he recurrido a otros códigos de internet, buscando aquellas variables que sirviesen para aumentar el score del modelo. En resumidas cuentas, lo que hacemos:
- Las variables con muchas categorías, como lo son
ward,funderyscheme_nametienen elementos, como paréntesis barras o signos de puntuación que hace que palabras iguales se cataloguen como distintas categorías. Por tanto, en primer lugar limpiamos el texto, y posteriormente hacemos lumping de igual que hemos hecho en clase. fe_popnum_ward: cantidad de poblacion por bomba de aguafe_dist: unifica las variables longitud y latitud (que son muy importantes para el modelo) en una nueva variable que codifique la distancia con respecto a un origen de coordenadas.fe_cant_agua: describe la cantidad de agua por población.fe_age: antigüedad de bomba de agua.fe_days: similar afe_agepero con respecto adate_recorded.fe_drdiff:año en el que se registró frente al que se construyó.
# cleanning some garbage present in ward and funder to reduce categories
clean_text <- function(text) {
stri_trans_tolower(
stri_replace_all_regex(
text,
pattern = "[ +\\p{Punct}()]",
replacement = "")
)
}
traintestimp <- traintestimp %>%
#population por bomb number * ward
#mutate(fe_popnum_ward = (sum(population) / n.()), .by = ward ) %>%
mutate(ward = clean_text(ward)) %>%
mutate(funder = clean_text(funder)) %>%
mutate(scheme_name = clean_text(scheme_name)) %>%
mutate(ward = as.character(fct_lump(ward, n = 200))) %>%
mutate(funder = as.character(fct_lump(funder, n = 200))) %>%
mutate(scheme_name = as.character(fct_lump(scheme_name, n = 200))) %>%
#geographic distance
mutate(fe_dist = distGeo(as.matrix(cbind(longitude, latitude)), c(0, 0))) %>%
# amount of water per population
mutate(fe_cant_agua = ifelse(population == 0, 0, round(amount_tsh / population, 3))) %>%
# age of the water pump
mutate(fe_age = 2014 - construction_year) %>%
mutate(fe_days = as.numeric(as.Date("2014-01-01") - as.Date(TrainTest$date_recorded))) %>%
# get the year and the month
mutate(fe_dryear = year(TrainTest$date_recorded)) %>%
mutate(fe_drmonth = month(TrainTest$date_recorded)) %>%
# difference between recorded date and construction date
mutate(fe_drdiff = fe_dryear-construction_year) %>%
as.data.table()Lo siguiente que hacemos es crear nuevas variables que expresen la frecuencia de las categorías para el caso de las variables categóricas.
cat_cols <- names(traintestimp %>% select(where(is.character)))
cat_cols## [1] "funder" "basin" "region"
## [4] "lga" "ward" "scheme_management"
## [7] "scheme_name" "extraction_type" "extraction_type_group"
## [10] "extraction_type_class" "management" "management_group"
## [13] "payment" "water_quality" "quality_group"
## [16] "quantity" "source" "source_type"
## [19] "source_class" "waterpoint_type" "waterpoint_type_group"
for (cat_col in cat_cols) {
# calculate for every value and then convert it to numeric
traintestimp[, paste0("fe_", cat_col) := as.numeric(.N), by = cat_col] }
# delete old cat cols
for (cat_col in cat_cols) {
traintestimp[, paste(cat_col) := NULL] }Con las variables de menos de 20 categorías, he probado a hacer dummies y crear modelos, pero el score de los mismos empeora con respecto a aquellos que no tienen las variables dummies.
datos_dummy <- dummy_cols(dummies) %>% select(where(is.numeric))
traintestimp <- cbind(traintestimp, datos_dummy)numeric_vars <- names(traintestimp %>% select(where(is.numeric)))
length(numeric_vars)## [1] 42
Como sabemos, algunos modelos son sensibles al escalado de los datos, por lo que una buena idea podría ser la de escalar los datos previo a ajustarlos a los distintos modelos.
# Función para estandarizar las variables
standardize <- function(x) {
(x - mean(x)) / sd(x) }
# Escalar las variables numéricas
traintestimp[, (numeric_vars) := lapply(.SD, standardize), .SDcols = numeric_vars]Sin embargo, esto no mejora en absoluto el rendimiento del modelo, por lo que optamos por prescindir de ello.
Hemos probado también a prescindir de ciertas variables que tenían poca importancia en los modelos de tipo random forest, de forma que de más importancia absoluta a otras variables que puedan tener más relevancia real y el score del modelo aumente, pero el hecho de eliminarlas no mejora el score real, por lo que optamos por dejarlas.
traintestimp %<>% select(-num_private) %>%
select(-NA_permit) %>%
select(-NA_public_meeting) %>%
select(-fe_dryear)Con todo esto, procedemos a dividir de nuevo el conjunto en
train y test para posteriormente crear y ajustar los
modelos. Para ello, creamos la columna idx.
# create idx column to separate
traintestimp %<>%
mutate( idx = 1:nrow(traintestimp)) %>%
as.data.table()
# separate train and test using idx.
xtrain <- traintestimp %>% as.data.table() %>%
filter( idx <= nrow(train)) %>%
select(-idx) %>%
#mutate( status_group = trainlab$status_group) %>% # we add the target
as.data.table()
xtest <- traintestimp %>% as.data.table() %>%
filter( idx > nrow(train)) %>%
select(-idx) %>%
as.data.table()Pasamos la variable objetivo vector_status_group a
numérico y posteriormente a factor, para poder realizar modelos como el
XGBoost o el lightGBM. Definimos la fórmula con la que
haremos los modelos.
vector_status_group <- trainlab$status_group
vector_status_group <- ifelse(vector_status_group == "functional", 0,
ifelse(vector_status_group == "functional needs repair", 1,
2))
vector_status_group <- as.factor(vector_status_group)
x_train <- xtrain %>% mutate(vector_status_group=vector_status_group) %>% as.data.frame()
x_test <- as.data.frame(xtest)
formula <- as.formula(vector_status_group ~ .)Podemos ver la clara desproporción de la variable objetivo, que tiene
muy pocas muestras de la clase 1. Ello hace que el modelo,
que ha sido entrenado con dichos datos, no sea capaz de captar la
información y, por tanto, a la hora de predecir tienda a fallar en dicha
clase. Por ello, es posible que la implementación de
SMOTEpara balancear las clases mejore el rendimiento del
modelo.
table(vector_status_group)## vector_status_group
## 0 1 2
## 32259 4317 22824
# completely unbalanced in class 1# Apply SMOTE
set.seed(0)
data <- x_train %>% mutate(vector_status_group = as.factor(vector_status_group))
# to increase the number of samples of the minoritary class
smote_data <- SMOTE(vector_status_group ~ ., data)#, perc.over = 100,perc.under = 500)
# Verify if smote worked correctly
table(smote_data$vector_status_group)##
## 0 1 2
## 10152 12951 7116
Pero la realidad es que el hecho de aplicar SMOTE no
solo no mejora el rendimiento si no que lo empeora, por lo que
prescindimos de ello.
3. Modelization
3.1. Model Selection
Una vez preprocesados los datos, realizamos un estudio mediante cross-validation para decidir cuál es el mejor modelo con el que ajustar a nuestros datos.
# Definir los parámetros
num_trees <- 100 # Número de árboles en el bosque
mtry <- 3 # Número de variables predictoras seleccionadas en cada división
n_folds <- 5 # Número de folds en la validación cruzada
# Crear un contenedor para almacenar los resultados de cada fold
accuracy_results <- numeric(n_folds)
# Realizar la validación cruzada
set.seed(42) # Fijar la semilla para reproducibilidad
#x_train <- xtrain %>% mutate(vector_status_group=vector_status_group) %>% as.data.frame()
# Obtener los índices de los folds
folds <- cut(seq(1, nrow(x_train)), breaks = n_folds, labels = FALSE)
rf_accuracy <- numeric(n_folds)
logistic_accuracy <- numeric(n_folds)
xgb_accuracy <- numeric(n_folds)
lgb_accuracy <- numeric(n_folds)
svm_accuracy <- numeric(n_folds)
# Iterar sobre cada fold
for (i in 1:n_folds) {
# Dividir los datos en conjuntos de entrenamiento y prueba según el fold actual
train_data <- x_train[folds != i, ]
test_data <- x_train[folds == i, ]
vector_train <- train_data %>% select(vector_status_group) %>% pull()
train_matrix <- train_data %>% select(-vector_status_group) %>% as.matrix()
vector_test <- test_data %>% select(vector_status_group) %>% pull()
test_matrix <- test_data %>% select(-vector_status_group) %>% as.matrix()
params <- list(
objective = "multiclass",
num_class = 3,
boosting_type = "gbdt",
metric = "multi_logloss")
# learning_rate = 0.02,
# max_depth = 15,
# num_threads = 4)
# Convertir los datos en formato xgb.DMatrix
lgb_train <- lgb.Dataset(data = train_matrix, label = as.numeric(vector_train) - 1)
lgb_test <- lgb.Dataset(data = test_matrix, label = as.numeric(vector_test) - 1)
# Ajustar el modelo de lightgbm en el conjunto de entrenamiento
capture.output({ lgbm <- lgb.train(data = lgb_train, nrounds = 600, valids = list(test = lgb_test), early_stopping_rounds = 10, params = params) }, file = NULL) #params = params,)
pred_lgb <- predict(lgbm, test_matrix)
pred_lgb_class <- max.col(matrix(pred_lgb, ncol = 3, byrow = TRUE)) - 1
lgb_accuracy[i] <- sum(pred_lgb_class == as.numeric(vector_test) - 1) / length(pred_lgb_class)
# Ajustar el modelo de Random Forest en el conjunto de entrenamiento
rf <- ranger(formula, importance = 'impurity', data = train_data, num.trees = num_trees, mtry = mtry, verbose = FALSE, seed = 42, probability = TRUE)
predrf <- predict(rf, data = test_data)$predictions
rf_accuracy[i] <- sum(max.col(predrf) == as.numeric(test_data$vector_status_group)) / length(test_data$vector_status_group)
# Ajustar el modelo de regresión logística multinomial en el conjunto de entrenamiento
glm <- suppressMessages(multinom(formula, data = train_data, MaxNWts = 10000, MaxIter = 100))
predglm <- predict(glm, newdata = test_data, type = "class")
logistic_accuracy[i] <- sum(predglm == test_data$vector_status_group) / length(test_data$vector_status_group)
# Ajustar el modelo de xgboost en el conjunto de entrenamiento
params <- list(
objective = "multi:softmax",
num_class = 3,
eval_metric = "mlogloss")
# Entrenar el modelo XGBoost
xgb_train <- xgb.DMatrix(data = train_matrix, label = as.numeric(vector_train) - 1)
xgb_test <- xgb.DMatrix(data = test_matrix, label = as.numeric(vector_test) - 1)
xgb_model <- xgb.train(params = params, data = xgb_train, nrounds = 600, nthread = 4)
predxgb <- predict(xgb_model, xgb_test)
xgb_accuracy[i] <- sum(predxgb == as.numeric(vector_test) - 1) / length(predxgb)
# Ajustar el modelo de svm en el conjunto de entrenamiento
train_data_svm <- data.frame(train_matrix, y = as.factor(as.numeric(vector_train) - 1))
test_data_svm <- data.frame(test_matrix, y = as.factor(as.numeric(vector_test) - 1))
svm_model <- svm(y ~ ., data = train_data_svm, method = "svmRadial", kernel = "radial", cost = 10, scale = TRUE)
pred_svm <- predict(svm_model, test_data_svm)
svm_accuracy[i] <- sum(pred_svm == as.numeric(vector_test) - 1) / length(pred_svm)
}3.2. Cross-Validation Accuracy per Model
Sin muchas dudas, el mejor modelo parece ser el XGBoost, ya que
es el que menor desviación tiene entre folds con el
cross-validation y es el que mayor valor medio de
accuracy tiene.
Aún siendo el XGBoost el mejor modelo, tuneamos además un random forest y creamos predicciones para ver cómo se desenvuelve y cuáles son las variables más importantes.
3.2. Random Forest Parameter Tunning
# Define the parameter grid to explore
grid <- expand.grid(
mtry = c(6,7,8), # Values for the number of predictor variables considered at each split
min.node.size = c(4, 5, 6), # Values for the minimum node size
splitrule = c("gini") # Values for the splitting rule
)
# Define the control function for hyperparameter search
ctrl <- trainControl(
method = "cv", # Use cross-validation
number = 5, # Number of cross-validation folds
search = "grid" # Perform grid search on the defined parameter grid
)
# Perform the hyperparameter search and train the model
rf_model <- train(
formula,
data = x_train,
method = "ranger",
trControl = ctrl,
tuneGrid = grid
)
# Show the best selected parameters
print(rf_model$bestTune)una vez encontrados los mejores parámetros, creamos el modelo y vemos el efecto de los parámetros en el score real de drivendata
3.2.1. num_trees effect in drivendata accuracy
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
3.2.2. mtry effect in drivendata accuracy
Con todo ello creamos el modelo y realizamos las predicciones.
# Train the final model with the best selected parameters
final_rf <- ranger(
formula,
importance = 'impurity',
data = x_train,
num.trees = 500, # Specify the number of trees directly here
mtry = 7,
min.node.size = 5,
splitrule = 'gini',
verbose = FALSE,
seed = 42,
probability = FALSE)
# predict test values
predi_rf <- predict(final_rf, x_test)
# change predicction from numeric to the categorical values
predi_rf <- ifelse(predi_rf$predictions == 0, "functional", ifelse(predi_rf$predictions == 1,"functional needs repair","non functional"))
# submit predictions
predi_rf <- data.table(id = test$id, status_group = predi_rf)
fwrite(predi_rf, file = "./submissions/subrf")Obtenemos una puntuación máxima de 0.8241 con el tuneado
del random forest.
3.2.3. Variable Importance
Por curiosidad, y en busca de obtener posible información, graficamos la importancia de las variables del modelo random forest, el cual es comúnmente usado como método de selección de variables importantes.
3.3. XGBoost Parameter Tunning
De igual forma, creamos y tuneamos un modelo XGBoost. En
primer lugar, convertimos los datos para poder trabajar con ellos.
train_matrix <- x_train %>% select(-vector_status_group) %>% as.matrix()
test_matrix <- xtest %>% as.matrix()
# Convert the data to DMatrix format
dtrain <- xgb.DMatrix(data = train_matrix, label = as.numeric(vector_status_group) - 1)
dtest <- xgb.DMatrix(data = test_matrix)hacemos el tuneado utilizando validación cruzada, de igual manera que en el random forest.
pb <- progress_bar$new(total = nrow(params_grid), format = "[:bar] :percent ETA: :eta")
# Define the hyperparameters for the grid search
params_grid <- expand.grid(
max_depth = c(6, 8, 10),
eta = c(0.01, 0.1),
gamma = c(0, 1),
colsample_bytree = c(0.5, 0.75, 1),
min_child_weight = c(1, 5),
subsample = c(0.5, 0.8),
objective = "multi:softprob",
num_class = length(unique(vector_status_group)) )
best_acc <- 0
best_params <- list()
pb <- progress_bar$new(total = nrow(params_grid), format = "[:bar] :percent ETA: :eta")
for (i in 1:nrow(params_grid)) {
params <- list(
max_depth = params_grid$max_depth[i],
eta = params_grid$eta[i],
gamma = params_grid$gamma[i],
colsample_bytree = params_grid$colsample_bytree[i],
min_child_weight = params_grid$min_child_weight[i],
subsample = params_grid$subsample[i],
objective = "multi:softprob",
num_class = length(unique(vector_status_group))
)
# Perform grid search for hyperparameter tuning
xgb_model <- xgb.cv(
params = params,
data = dtrain,
nfold = 5,
nrounds = 200,
early_stopping_rounds = 20,
verbose = 0
)
if (1 - tail(xgb_model$evaluation_log$test_mlogloss_mean, 1) > best_acc) {
best_acc <- 1 - tail(xgb_model$evaluation_log$test_mlogloss_mean, 1)
best_params <- params
}
# progress bar to show the %
pb$tick()
}
print(best_params)de forma similar a lo que hemos hecho antes, creamos el modelo y
vemos el efecto del parámetro rounds en el score
real de drivendata, ya que es el que más influye en el
overfitting y el que más influye en el tiempo de ejecución, por
lo que no lo he usado como parámetro de entrada en el
tunning
3.3.1. rounds effect in drivendata accuracy
una vez tenemos las parámetros ajustados, creamos el modelo y las predicciones
params <- list(
max_depth = 15,
eta = 0.02,
gamma = 0,
colsample_bytree = 0.3,
colsample_bylevel = 0.9,
colsample_bynode = 0.7,
min_child_weight = 1,
subsample = 0.8,
objective = "multi:softmax",
num_class = 3)
final_xgb <- xgboost(
params = params,
data = dtrain,
nrounds = 600,
verbose = 0)
# Make predictions on the test dataset
predi_xgb <- predict(final_xgb,dtest)
predi_xgb <- ifelse(predi_xgb == 0, "functional", ifelse(predi_xgb == 1, "functional needs repair", "non functional"))
# Submit predictions
predi_xgb <- data.table(id = test$id, status_group = predi_xgb)
fwrite(predi_xgb, file = "./submissions/subxgb")obtenemos la que hasta el momento es la mejor score.
3.3.2. Variable Importance
3.3. Stacking
Por último, recurrimos a una de las técnicas más utilizadas en los concursos, el stacking. Con nuestros modelos tuneados y las predicciones hechas, tratamos de mejorar el resultado combinando el random forest con el XGBoost.
split_index <- createDataPartition(x_train$vector_status_group, p = 0.7, list = FALSE)
x_train_base <- x_train[split_index, ]
x_train_meta <- x_train[-split_index, ]
# Convertir x_train_base y x_train_meta a matrices para XGBoost
y_train_base <- vector_status_group[split_index]
y_train_meta <- vector_status_group[-split_index]
x_train_base_matrix <- x_train_base %>% select(-vector_status_group) %>% as.matrix()
x_train_meta_matrix <- x_train_meta %>% select(-vector_status_group) %>% as.matrix()# Ranger model
ranger_model <- ranger(formula,
data = x_train_base, num.trees = 500,
mtry = 7, importance = 'impurity',
min.node.size = 5, splitrule = 'gini',
verbose = FALSE, seed = 42, probability = FALSE)
# Predecir con el modelo Ranger
ranger_predictions_meta <- predict(ranger_model, x_train_meta, type = "response")$predictions# Xgboost model
params <- list(max_depth = 15, eta = 0.02,
gamma = 0, colsample_bytree = 0.3,
colsample_bylevel = 0.9, colsample_bynode = 0.7,
min_child_weight = 1, subsample = 0.8,
objective = "multi:softmax", num_class = 3)
# Convert the data to DMatrix format
dtrain_base <- xgb.DMatrix(data = x_train_base_matrix, label = as.numeric(y_train_base) - 1)
# train model
xgb_model <- xgb.train(params = params, data = dtrain_base, nrounds = 600)
# Predecir con el modelo XGBoost
xgb_predictions_meta <- predict(xgb_model, x_train_meta_matrix)Ahora, creamos el modelo con el que hacemos stacking. Elegimos la
regresión lineal como segunda opción ya que el randomForest
no ha logrado mejorar el score.
# Crear un data.frame con las predicciones de ambos modelos
combined_predictions_meta <- data.frame(ranger = ranger_predictions_meta, xgboost = xgb_predictions_meta)
# Entrenar el modelo meta
#meta_model <- randomForest(y ~ ., data = cbind(y = x_train_meta$vector_status_group, combined_predictions_meta), ntree = 500)
library(nnet)
meta_model <- multinom(y ~ ., data = cbind(y = x_train_meta$vector_status_group, combined_predictions_meta))# Predecir con el modelo Ranger
ranger_predictions_test <- predict(ranger_model, x_test, type = "response")$predictions
# Predecir con el modelo XGBoost
test_matrix <- xtest %>% as.matrix()
xgb_predictions_test <- predict(xgb_model, test_matrix)# Crear un data.frame con las predicciones de ambos modelos
combined_predictions_test <- data.frame(ranger = ranger_predictions_test, xgboost = xgb_predictions_test)
# Stacking de los modelos Ranger y XGBoost
stacked_predictions <- predict(meta_model, combined_predictions_test)
# Make predictions on the test dataset
stacked_predictions <- ifelse(stacked_predictions == 0, "functional", ifelse(stacked_predictions == 1, "functional needs repair", "non functional"))
# Submit predictions
stacked_predictions <- data.table(id = test$id, status_group = stacked_predictions)
fwrite(stacked_predictions, file = "./submissions/stacking")Hacer stacking empeora el score del modelo XGBoost,
obteniendo un score de 0.816.
Por tanto, nos quedamos con el que hasta el momento es el mejor
modelo, el XGBoost.
driven data score